Przed oddaniem zadania upewnij się, że wszystko działa poprawnie. Uruchom ponownie kernel (z paska menu: Kernel$\rightarrow$Restart) a następnie wykonaj wszystkie komórki (z paska menu: Cell$\rightarrow$Run All).
Upewnij się, że wypełniłeś wszystkie pola TU WPISZ KOD lub TU WPISZ ODPOWIEDŹ, oraz
że podałeś swoje imię i nazwisko poniżej:
NAME = "Piotr Durniat"
W tym notebooku zawarte są zadania wprowadzające do rachunku prawdopodobieństwa i statystyki w środowisku python oraz programowania probabilistycznego.
Rachunek prawdopodobieństwa, czyli dział matematyki zajmujący się niepewnością i losowością, stanowi podstawę rozpoznawania wzorców. Umożliwia on ujęcie ilościowe niepewności wynikającej np. z szumów występujących w danych lub też ograniczonej ich ilości, dzięki czemu możemy wykonywać predykcje na podstawie danych, które są niekompletne lub niejednoznaczne.
W zadaniu będziemy posługiwać się przykładem. Wyobraźmy sobie dwa pudełka: jedno czerwone i jedno niebieskie. Czerwone pudełko zawiera w sobie 4 jabłka i 6 pomarańczy, natomiast niebieskie zawiera 7 jabłek i 3 pomarańcze. Załóżmy, że:
W tym przykładzie, wybór pudełka jest zmienną losową; oznaczmy ją jako $B$. Zmienna ta może przyjmować dwie wartości: $B=r$ (czerwone pudełko) albo $B=b$ (niebieskie pudełko). Analogicznie, wybór owocu jest również zmienną losową; oznaczmy ją jako $F$ - wówczas $F=a$ oznacza wybór jabłka, a $F=o$ oznacza wybór pomarańczy.
Prawdopodobieństwo wyboru czerwonego pudełka wynosi $\Pr(B=r) = 0.3$, natomiast prawdopodobieństwo wyboru niebieskiego pudełka to $\Pr(B=b) = 0.7$. Są to zdarzenia rozłączne oraz obejmują wszystkie możliwe wyniki, zatem suma prawdopodobieństw ich wystąpienia wynosi $\Pr(B=r) + \Pr(B=b) = 1$.
Załóżmy, że powtarzamy powyższy schemat czterokrotnie w celu wylosowania czterech owoców (kolejność nie jest istotna). Podaj przestrzeń zdarzeń elementarnych dla doświadczenia losowego. Zapisz ją w formie listy krotek odpowiednich znaków (np. [("a", "b"), ("c", "d")] i przypisz do zmiennej space_1.
# "a" - apple
# "o" - orange
# kolejność nie istotna - liczy się tylko liczba jabłek i pomarańczy w krotce
space_1 = [
("a", "a", "a", "a"),
("a", "a", "a", "o"),
("a", "a", "o", "o"),
("a", "o", "o", "o"),
("o", "o", "o", "o"),
]
sample_space_1_test
Score: 0.25 / 0.25 (Top)
assert type(space_1) == list
for sample in space_1:
assert type(sample) == tuple
assert len(sample) == 4
for item in sample:
assert type(item) == str
assert len(item) == 1
assert item in {"a", "o", "b", "r"}
### POCZĄTEK UKRYTYCH
sorted_space_1 = [sorted(sample) for sample in space_1]
assert sorted(('a', 'a', 'a', 'a')) in sorted_space_1
assert sorted(('a', 'a', 'a', 'o')) in sorted_space_1
assert sorted(('a', 'a', 'o', 'o')) in sorted_space_1
assert sorted(('a', 'o', 'o', 'o')) in sorted_space_1
assert sorted(('o', 'o', 'o', 'o')) in sorted_space_1
### KONIEC UKRYTYCH
Drugie z doświadczeń polega na wylosowaniu pary ${B, F}$: najpierw losujemy pudełko, a następnie losujemy z niego owoc. Podaj przestrzeń zdarzeń dla tak zdefiniowanego doświadczenia losowego. Zapisz je w postaci listy krotek (np. [("a", "b"), ("c", "d")] i przypisz do zmiennej space_2.
# B (Box):
# - "b" - blue
# - "r" - red
# F (Fruit):
# - "a" - apple
# - "o" - orange
space_2 = [
("b", "a"),
("b", "o"),
("r", "a"),
("r", "o"),
]
sample_space_2_test
Score: 0.25 / 0.25 (Top)
assert type(space_2) == list
for sample in space_2:
assert type(sample) == tuple
assert len(sample) == 2
assert sample[0] in {"r", "b"}
assert sample[1] in {"a", "o"}
assert all(type(item) == str for item in sample)
### POCZĄTEK UKRYTYCH
assert ('r', 'a') in space_2
assert ('r', 'o') in space_2
assert ('b', 'a') in space_2
assert ('b', 'o') in space_2
### KONIEC UKRYTYCH
Uogólniając powyższy przykład, rozważmy dwie zmienne losowe $X$ i $Y$. Zmienna $X$ może przyjąć dowolną wartość $x \in \{x_1, x_2, \ldots, x_M\}$, a zmienna $Y$ może przyjąć wartości $y \in \{y_i, y_2, \ldots, y_L\}$. Rozważmy, że $N$ razy próbkujemy wartości obu zmiennych $X$ i $Y$; liczba prób, w których $X=x_i$ i $Y=y_j$ wynosi $n_{ij}$. Ponadto, załózmy, że liczba przypadków, gdzie $X$ przyjmuje wartość $x_i$ (bez względu na wartość $Y$) to $c_i$; analogicznie liczba prób, gdy $Y=y_j$ oznaczamy przez $r_j$.
Prawdopodobieństwo łączne, tzn. prawdopodobieństwo, że zmienna $X$ przyjmie wartość $x_i$ oraz że $Y$ przyjmie wartość $y_j$ definiuje się jako:
$$\Pr(X=x_i, Y=y_j) = \frac{n_{ij}}{N}\tag{1}$$Implicite rozważamy tu granicę przy $N\rightarrow\infty$. Analogicznie, prawdopodobieństwo brzegowe, tzn. prawdopodobieństwo, że $X$ przyjmuje wartość $x_i$ bez względu na wartość $Y$ można zapisać jako:
$$\Pr(X=x_i) = \frac{c_i}{N}\tag{2}$$W doświadczeniu losowym z pudełkami i owocami, rozpisz liczby wystąpień każdego z przypadków w tablicy dwudzielczej, które dadzą określone wcześniej prawdopodobieństwa, przekształcając powyższe wzory i przypisując odpowiednie wartości do zmiennych zadeklarowanych w tablicy. Tablica dwudzielcza to tabela przedstawiająca łączny rozkład dwóch zmiennych: każda komórka zawiera liczbę wystąpień kombinacji wartości zmiennych. Załóż, że wykonano $N=100$ prób.
| B / F | B=r | B=b | |
|---|---|---|---|
| F=a | n_ra |
n_ba |
$$\sum_B=r_a$$ |
| F=o | n_ro |
n_bo |
$$\sum_B=r_o$$ |
| $$\sum_F=c_r$$ | $$\sum_F=c_b$$ | $$\sum_B\sum_F = N$$ |
p_r = 0.3
p_b = 0.7
# Probability of picking an apple if the box is red
p_a_if_r = 0.4
p_o_if_r = 0.6
p_a_if_b = 0.7
p_o_if_b = 0.3
n = 100
n_ra = p_r * p_a_if_r * n
n_ro = p_r * p_o_if_r * n
n_ba = p_b * p_a_if_b * n
n_bo = p_b * p_o_if_b * n
crosstab_tests
Score: 0.25 / 0.25 (Top)
assert n_ra + n_ro + n_ba + n_bo == 100
### POCZĄTEK UKRYTYCH
cr = n_ra + n_ro
cb = n_ba + n_bo
ra = n_ra + n_ba
ro = n_ro + n_bo
assert cr == 30
assert cb == 70
assert n_ra == 12
assert n_ro == 18
assert n_ba == 49
assert n_bo == 21
### KONIEC UKRYTYCH
Marginalizacja polega na przejściu z rozkładu łącznego na rozkład brzegowy. Patrząc na tablicę dwudzielną, polega ona na zsumowaniu liczby wystąpień w wierszach albo kolumnach i podzieleniu przez całkowitą liczbę próbek (reguła brzegowa).
$$\Pr(X=x_i) = \sum_{j=1}^L \Pr(X=x_i, Y=y_j)\tag{3}$$Jeżeli rozważymy przykłady, dla których $X=x_i$, wówczas część spośród nich, dla których $Y=y_j$ określa się mianem prawdopodbieństwa warunkowego:
$$\Pr(Y=y_j \mid X=x_i) = \frac{n_{ij}}{c_i}\tag{4}$$Stąd też wynika zależność - reguła łańcuchowa:
$$\Pr(X=x_i, Y=y_j) = \frac{n_{ij}}{N} = \frac{n_{ij}}{c_i} \frac{c_i}{N} = \Pr(Y=y_j \mid X=x_i)\Pr(X=x_i)\tag{5}$$Oblicz wartości prawdopodobieństwa warunkowego dla wcześniejszego doświadczenia losowego:
| prawdopodobieństwo | zmienna |
|---|---|
| $$\Pr(F=a \mid B=r)$$ | p_fa_br |
| $$\Pr(F=o \mid B=r)$$ | p_fo_br |
| $$\Pr(F=a \mid B=b)$$ | p_fa_bb |
| $$\Pr(F=o \mid B=b)$$ | p_fo_bb |
| $$\Pr(B=r \mid F=a)$$ | p_br_fa |
| $$\Pr(B=b \mid F=a)$$ | p_bb_fa |
| $$\Pr(B=r \mid F=o)$$ | p_br_fo |
| $$\Pr(B=b \mid F=o)$$ | p_bb_fo |
c_r = n_ra + n_ro
c_b = n_ba + n_bo
r_a = n_ra + n_ba
r_o = n_ro + n_bo
p_fa_br = n_ra / c_r
p_fo_br = n_ro / c_r
p_fa_bb = n_ba / c_b
p_fo_bb = n_bo / c_b
p_br_fa = n_ra / r_a
p_bb_fa = n_ba / r_a
p_br_fo = n_ro / r_o
p_bb_fo = n_bo / r_o
conditional_probs_test
Score: 0.25 / 0.25 (Top)
# testy ukryte
### POCZĄTEK UKRYTYCH
from nose.tools import assert_almost_equal
assert_almost_equal(p_fa_br, 0.400, places=3)
assert_almost_equal(p_fo_br, 0.600, places=3)
assert_almost_equal(p_fa_bb, 0.700, places=3)
assert_almost_equal(p_fo_bb, 0.300, places=3)
assert_almost_equal(p_br_fa, 0.197, places=3)
assert_almost_equal(p_bb_fa, 0.803, places=3)
assert_almost_equal(p_br_fo, 0.462, places=3)
assert_almost_equal(p_bb_fo, 0.538, places=3)
### KONIEC UKRYTYCH
Przekształcając równanie reguły łańcuchowej $(5)$ oraz reguły brzegowej $(3)$ oraz wykorzystując właściwość, że $\Pr(X, Y) = \Pr(Y, X)$, uzyskujemy zależność pomiędzy prawdopodobieństwami warunkowymi:
$$\Pr(Y \mid X) = \frac{\Pr(X \mid Y) \Pr(Y)}{\Pr(X)}\tag{6}$$Zależność ta nazywana jest regułą Bayesa i pełni ona centralną rolę w probabilistycznym uczeniu maszynowym. Wykorzystując ponownie regułę brzegową $(3)$, mianownik równania $(6)$ można przedstawić w postaci:
$$\Pr(X) = \sum_Y \Pr(X \mid Y) \Pr(Y)\tag{7}$$Można go zatem potraktować jako "stałą normalizującą", gwarantującą, że prawdopodobieństwo warunkowe $\Pr(Y \mid X)$ sumuje się do $1$.
Test diagnostyczny dotyczący pewnej choroby ma następujące parametry:
Zapisz wynik w zmiennej p_y1_x1.
Z prawa całkowitego prawdopodobieństwa:
Z prawa Bayesa:
$$ P(Y=1\mid X=1) = $$$$ = \frac{P(X=1 \mid Y=1) P(Y=1)}{P(X=1)} = $$$$ = \frac{P(X=1 \mid Y=1) P(Y=1)} {\sum_Y P(X=1 \mid Y) P(Y)} = $$$$ = \frac{P(X=1 \mid Y=1) P(Y=1)} {P(X=1 \mid Y=1) P(Y=1) + P(X=1 \mid Y=0) P(Y=0)} = $$$$ = \frac{P(X=1 \mid Y=1) P(Y=1)} {P(X=1 \mid Y=1) P(Y=1) + (1 - P(X=0\mid Y=0))(1-P(Y=1))} = $$$$ = \frac{0.9 \cdot 0.014}{0.9 \cdot 0.014 + (1-0.03) \cdot (1-0.014)} \approx 0.299 $$p_x1_y1 = 0.9
p_y1 = 0.014
p_x0_y0 = 0.97
p_y1_x1 = p_x1_y1 * p_y1 / (p_x1_y1 * p_y1 + (1 - p_x0_y0) * (1 - p_y1))
print(p_y1_x1)
bayes_prob
Score: 0.25 / 0.25 (Top)
# testy ukryte
### POCZĄTEK UKRYTYCH
from nose.tools import assert_almost_equal
assert_almost_equal(p_y1_x1, 0.299, places=3)
### KONIEC UKRYTYCH
Oprócz opisywanych wcześniej zmiennych losowych dyskretnych, prawdopodobieństwa mogą dotyczyć również zmiennych ciągłych. Rozkład prawdopodobieństa zmiennych ciągłych opisany jest funkcją gęstości $p(x)$. Prawdopodobieństwo, że zmienna losowa $X$ przyjmuje wartości z przedziału $[a, b]$ można wyznaczyć przy pomocy równania:
$$\Pr(a \le X \le b) = \int_a^b p(x) dx\tag{8}$$Funkcja gęstości spełnia dwa kryteria:
Dystrybuantą, czyli funkcją skumulowanej gęstości prawdopodobieństwa nazywamy prawdopodobieństwo, że $x$ leży w przedziale $(-\infty, z)$:
$$F(z) = \int_{-\infty}^z p(x) dx\tag{9}$$Dla zmiennych dyskretnych funkcja $p(x)$ jest nazywana funkcją masy prawdopodobieństwa.
Załóżmy, że $x$ i $y$ są dwoma zmiennymi ciągłymi. Funkcja gęstości spełnia reguły:
Podstawową własnością rozkładu jest wartość oczekiwana, określająca spodziewany wynik doświadczenia losowego. Dla rozkładu dyskretnego dany jest wzorem:
$$\mathbb{E}[X] = \sum_{x\in\mathcal{X}} x \cdot p(x)\tag{10}$$natomiast dla rozkładu ciągłego definiuje się:
$$\mathbb{E}[X] = \int_{x\in\mathcal{X}} x \cdot p(x) dx\tag{11}$$Wariancja to własność określająca jak bardzo wartości $x$ odchylają się od średniej (rozpiętość rozkładu)
$$Var[X] = \mathbb{E}[(X - \mathbb{E}[X])^2] = \mathbb{E}[X^2] - (\mathbb{E}[X])^2\tag{12}$$Odchylenie standardowe definiowane jest jako pierwiastek z wariancji:
$$std[X] = \sqrt{Var[X]}\tag{13}$$Dla dwóch zmiennych losowych $X$ i $Y$, kowariancja określa stopień, w jakim zależą od siebie liniowo. Dana jest ona wzorem:
$$cov[X, Y] = \mathbb{E}[(X - \mathbb{E}[X])(Y - \mathbb{E}[Y])] = \mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y] \in [0, \infty)\tag{14}$$Korelacja jest znormalizowana postaciÄ… kowariancji:
$$corr[X, Y] = \frac{cov[X, Y]}{\sqrt{Var[X]Var[Y]}} \in [-1, 1]\tag{15}$$Dla liniowo zależnych zmiennych losowych $X$ i $Y$ $|corr[X, Y]| = 1$. Gdy zmienne są od siebie niezależne, tzn. $\Pr(X,Y) = \Pr(X)\Pr(Y)$, wtedy $corr[X, Y] = 0$.
Pyro to probabilistyczny język programowania udostępniony przez Uber AI Labs. Zbudowany jest on na frameworku PyTorch, łącząc modelowanie modeli głębokich z modelowaniem bayesowskim.
Probabilistyczny język programowania to język zaprojektowany do budowy oraz wnioskowania w modelach probabilistycznych. Zagadnienie to będzie omawiane w ramach kolejnych laboratoriów - tutaj skupimy się na podstawowych prymitywach Pyro.
Podobnie jak PyTorch, Pyro pozwala na wykorzystanie dowolnych instrukcji Pythona, takich jak pętle, rekursja, funkcje wysokopoziomowe itp., dzięki czemu można w nim zareprezentować dowolny obliczalny rozkład prawdopodobieństwa.
Podstawowym elementem programów probabilistycznych są funkcje stochastyczne. W Pyro są to pythonowe callable (obiekty implementujące metodę __call__()) lub moduły PyTorcha nn.Module, zawierające deterministyczny kod oraz podstawowe funkcje stochastyczne, wywołujące generator liczb losowych. W ramach tego laboratorium rozważać będziemy te pierwsze.
Pyro wykorzystuje funkcje stochastyczne (rozkłady prawdopodobieństwa) będące wrapperem torch.distribution. Znajdują się one w module pyro.distributions (dokumentacja). Aby próbkować wartości z tych rozkładów można używać metod z biblioteki torch: .rsample() oraz .sample(), jednak w przypadku Pyro lepiej wykorzystać funkcję pyro.sample, będącą jedną z kluczowych elementów biblioteki. Jej użycie przedstawiono w kolejnej komórce.
import pyro.distributions as dist
distribution = dist.Normal(0, 1)
# te dwa wywołania są równoważne
torch_sampled = distribution.sample()
# pyro.sample pozwala na tworzenie programu Pyro
pyro_sampled = pyro.sample("my_sample", distribution)
# przykład użycia pyro.distributions dla rozkładu kategorycznego i normalnego
from IPython.display import display, Math
import pyro
import pyro.distributions as dist
import torch
# rozkład kategoryczny
probas = [0.2, 0.15, 0.33, 0.26, 0.06] # prawdopodobieństwa kategorii
categorical = dist.Categorical(torch.tensor(probas)) # rozkład kategoryczny
display(Math(f"x \sim Cat({probas})"))
x = pyro.sample( # próbkowanie z rozkładu
"category", # nazwa zmiennej przechowywanej przez pyro
categorical, # podstawowa funkcja stochastyczna
)
x_log_prob = categorical.log_prob(x) # prawdopodobieństwo zaobserwowania
# tej wartości przy próbkowaniu
# z rozkładu kategorycznego
display(Math(f"\Pr(x = {x}) = {torch.exp(x_log_prob).item():.3f})"))
# rozkład normalny
loc = 0.0 # wartość oczekiwana 0.
scale = 1.0 # odchylenie standardowe 1.
normal = dist.Normal(loc, scale) # rozkład normalny
display(Math(f"y \sim \mathcal{{N}}(\mu={loc}, \sigma={scale})"))
y = pyro.sample("normal", normal) # próbkowanie z rozkładu
y_log_prob = normal.log_prob(y) # prawdopodobieństwo zaobserwowania
# tej wartości przy próbkowaniu
# z rozkładu normalnego
display(Math(f"\Pr(y = {y:.3f}) = {torch.exp(y_log_prob).item():.3f})"))
Wybierz jeden ciągły rozkład prawdopodobieństwa spośród dostępnych w pyro.distributions i omówionych na wykładzie. Przygotuj animację, która pokaże jak z kolejnymi próbkowaniami rozkładu zmieniać się będzie rozkład wartości. Stwórz histogram wystąpień.
Podpowiedź: próbkowane wartości zbieraj w liście i generuj wykres dla każdego kroku. Skorzystaj z przygotowanej funkcji animowania.
## przykład tworzenia animacji z użyciem matplotlib.animation
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import HTML
from matplotlib import animation, rc
def animate():
# przygotuj wykres do rysowania animacji
fig, ax = plt.subplots()
ax.set_xlim(( 0, 2))
ax.set_ylim((-2, 2))
line, = ax.plot([], [], lw=2)
# funkcja animacji, wywoływana sekwencyjnie
def sin(i):
x = np.linspace(0, 2, 1000)
y = np.sin(2 * np.pi * (x - 0.01 * i))
line.set_data(x, y)
return (line,)
# funkcja inicjalizacji, rysująca tło każdej klatki
def init():
line.set_data([], [])
return (line,)
# wywołanie animacji; blit=True powoduje rysowanie tylko zmian
anim = animation.FuncAnimation(
fig, sin, init_func=init, frames=100, interval=20, blit=True
)
return anim.to_jshtml()
HTML(animate())
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from pyro.distributions import Exponential
# Inicjalizacja rozkładu wykładniczego
lambda_param = 0.7
exponential_dist = Exponential(rate=lambda_param)
def animate_distribution():
fig, ax = plt.subplots()
# Ustawienie zakresów osi
ax.set_xlim((0, 10))
ax.set_ylim((0, 1.5))
data = []
(hist,) = ax.plot([], [], lw=2)
# Funkcja rysujÄ…ca kolejne klatki animacji
def draw(i):
sample = exponential_dist.sample().item()
data.append(sample)
ax.cla()
ax.hist(data, bins=30, range=(0, 10), density=True)
ax.set_xlim((0, 10))
ax.set_ylim((0, 1.5))
return (hist,)
# Funkcja inicjalizujÄ…ca
def init():
ax.hist([], bins=30, range=(0, 10), density=True)
return (hist,)
# Wywołanie animacji
anim = animation.FuncAnimation(
fig, draw, init_func=init, frames=200, interval=100, blit=False
)
return anim.to_jshtml()
from IPython.display import HTML
HTML(animate_distribution())
Tak jak wspomniano wcześniej, w Pyro funkcje stochastyczne mogą zawierać zarówno odwołania do podstawowych funkcji stochastycznych pyro.distributions, jak i deterministyczny kod pythonowy, co pozwala na modelowanie różnych doświadczeń losowych. Przykładowo, w poniższej komórce znajduje się model doświadczenia losowego, przedstawionego na początku notatnika.
# doświadczenie losowe: wylosuj pudełko, a następnie z niego wylosuj owoc
def get_fruit(box_proba=0.45, red_box_proba=0.4, blue_box_proba=0.7):
boxes = ["r", "b"]
fruit = ["a", "o"]
box_idx = pyro.sample("box", dist.Bernoulli(box_proba))
if box_idx == 0:
fruit_idx = pyro.sample("fruit", dist.Bernoulli(red_box_proba))
else:
fruit_idx = pyro.sample("fruit", dist.Bernoulli(blue_box_proba))
return boxes[int(box_idx)], fruit[int(fruit_idx)]
for _ in range(5):
box, fruit = get_fruit()
print(f"Wylosowano '{fruit}' z pudełka '{box}'")
Dzięki temu możliwe jest stworzenie symulatora, który powtórzy doświadczenie dostatecznie dużą liczbę razy w celu numerycznego przybliżenia prawdopodobieństw, które wcześniej wyznaczyliśmy obliczając je ręcznie.
# powtórz doświadczenie dużą liczbę razy i aproksymuj prawdopodobieństwa
total = 10_000
n_boxes = {"r": 0, "b": 0}
n_fruit = {"a": 0, "o": 0}
for _ in range(total):
box, fruit = get_fruit()
n_boxes[box] += 1
n_fruit[fruit] += 1
display(Math(f"$$\Pr(B=r)\simeq{n_boxes['r'] / total:.3f}$$"))
display(Math(f"$$\Pr(B=b)\simeq{n_boxes['b'] / total:.3f}$$"))
display(Math(f"$$\Pr(F=a)\simeq{n_fruit['a'] / total:.3f}$$"))
display(Math(f"$$\Pr(F=o)\simeq{n_fruit['o'] / total:.3f}$$"))
print("itd.")
Rozważmy doświadczenie losowe: student aplikował na 9 kierunków. Szanse przyjęcia dla każdego z kierunków są równe i wynoszą 1%. Niestety, student nie został przyjęty na żaden z kierunków. Jakie jest prawdopodobieństwo takiej sytuacji? Wybierz odpowiedni rozkład prawdopodobieństwa i oblicz prawdopodobieństwo ręcznie, wpisując tok rozumowania do następnej komórki oraz przypisując wynik do zmiennej p_reject. Zamodeluj doświadczenie w formie symulatora uzupełniając funkcję attempt_join_university, która wykona doświadczenie losowe, oraz funkcję simulate_join_university, która wykona to doświadczenie zadaną liczbę razy i zwróci prawdopodobieństwo opisywanego zdarzenia.
Prawdopodobieństwo odrzucenia na wszystkich kierunkach:
$$ p_{reject} = (1-p_{admission})^9 = 0.99^9 = 0.91352 $$p_join_one = 0.01
p_reject_one = 1 - p_join_one
application_count = 9
p_reject = p_reject_one**application_count
admission_manual_tests
Score: 0.25 / 0.25 (Top)
display(Math(f"$$\Pr(reject)={p_reject:.5f}$$"))
# testy ukryte
### POCZĄTEK UKRYTYCH
from nose.tools import assert_almost_equal
assert_almost_equal(p_reject, 0.914, places=3)
### KONIEC UKRYTYCH
def attempt_joint_faculty(p_admission):
"""
Funkcja modelująca doświadczenie losowe dla jednego kierunku.
"""
return pyro.sample("application", dist.Bernoulli(p_admission))
def attempt_join_university(n_trials=9, p_admission=0.01):
"""
Funkcja stochastyczna modelująca doświadczenie losowe.
Zwraca 1 jeśli student został odrzucony na wszystkie kierunki.
"""
return all(attempt_joint_faculty(p_admission) == 0 for _ in range(n_trials))
def simulate_join_university(n_simulations=10_000):
"""Wykonaj doświadczenie losowe określoną liczbę.
Zwraca przybliżone prawdopodobieństwo nie przyjęcia na żaden z kierunków
przy zadanej liczbie prób.
"""
n_rejected = 0
for _ in range(n_simulations):
if attempt_join_university() == 1:
n_rejected += 1
return n_rejected / n_simulations
admission_code_tests
Score: 1.0 / 1.0 (Top)
display(Math(f"$$\Pr(reject)\simeq{simulate_join_university():.5f}$$"))
# testy ukryte
### POCZĄTEK UKRYTYCH
from nose.tools import assert_almost_equal
assert_almost_equal(simulate_join_university(100), 0.914, places=1)
assert_almost_equal(simulate_join_university(1_000), 0.914, places=1)
assert_almost_equal(simulate_join_university(10_000), 0.914, places=1)
### KONIEC UKRYTYCH
Przygotuj symulator dla zadania 1.4. Sprawdź, jak czułość i swoistość testu oraz prawdopodobieństwo choroby wpływają na osiągany wynik. Wykonaj wizualizację zależności prawdopodobieństwa bycia chorym przy pozytywnym wyniku testu od ww. parametrów.
import numpy as np
import matplotlib.pyplot as plt
import pyro.distributions as dist
def p_illness_given_positive(sensitivity, specificity, p_ilness):
# (P(X=1|Y=1)P(Y=1)} / (P(X=1|Y=1)P(Y=1)+(1-P(X=0|Y=0))(1-P(Y=1))) =
return (
sensitivity
* p_ilness
/ (sensitivity * p_ilness + (1 - specificity) * (1 - p_ilness))
)
def draw_chart():
sensitivity_range = np.linspace(0.1, 1.0, 10)
specificity_range = np.linspace(0.1, 1.0, 10)
disease_probability_range = np.linspace(0.01, 0.99, 10)
S, D, P = np.meshgrid(
sensitivity_range, specificity_range, disease_probability_range
)
PPV = p_illness_given_positive(S, D, P)
print(PPV.shape)
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection="3d")
ax.scatter(S.flatten(), D.flatten(), P.flatten(), c=PPV.flatten(), cmap="viridis")
ax.set_xlabel("Sensitivity")
ax.set_ylabel("Specificity")
ax.set_zlabel("Positive Predictive Value")
ax.set_title("PPV vs Sensitivity and Specificity")
def simulate_illness(sensitivity: float, specificity: float, p_ilness: float):
"""
Simulate the test result
return test result and actual illness status
"""
y = pyro.sample("test", dist.Bernoulli(p_ilness))
if y == 1:
return pyro.sample("ilness", dist.Bernoulli(sensitivity)).item(), y
else:
return pyro.sample("ilness", dist.Bernoulli(1 - specificity)).item(), y
def simulate_p_illness_given_positive(
sensitivity, specificity, p_illness, n_simulations=1_000
):
positive_tests = 0
true_positives = 0
for _ in range(n_simulations):
x, y = simulate_illness(sensitivity, specificity, p_illness)
if x == 1: # If the test result is positive
positive_tests += 1
if y == 1: # If the person is actually ill
true_positives += 1
return (true_positives / positive_tests) if positive_tests > 0 else 0
# Czułość P(X=1|Y=1)
sensitivity = 0.9
# Swoistość (P(X=0|Y=0)
specificity = 0.97
# P(Y=1)
p_ilness = 0.014
print("Calculated P: ", p_illness_given_positive(sensitivity, specificity, p_ilness))
print(
"Simulated P: ",
simulate_p_illness_given_positive(sensitivity, specificity, p_ilness),
)
draw_chart()